Momentum Transport in Granular Flows 



Gregg LoisW, Anael Lemaitre^ 2 -', and Jean M. Carlson^ 1 ) 
^ Department of Physics, University of California, Santa Barbara, California 93106, U.S.A. and 
Institut Navier- LMSGC, 2 allee Kepler, 77420 Champs-sur-Marne, France 
(Dated: February 5, 2008) 

We investigate the error induced by only considering binary collisions in the momentum transport 
of hard-sphere granular materials, as is done in kinetic theories. In this process, we first present a 
general microscopic derivation of the momentum transport equation and compare it to the kinetic 
theory derivation, which relies on the binary collision assumption. These two derivations yield 
different microscopic expressions for the stress tensor, which we compare using simulations. This 
provides a quantitative bound on the regime where binary collisions dominate momentum transport 
and reveals that most realistic granular flows occur in the region of phase space where the binary 
collision assumption does not apply. 

Granular materials have emerged as an ideal non-equilibrium system, offering the opportunity to test statistical 
physics in settings where the underlying assumptions can be closely examined 0, 0, A central objective for 
granular flows is the formation of a hydrodynamic description, which is essential for both theories and applications. 
To achieve this goal, there has been considerable interest in describing granular flows using extensions of the kinetic 
theory of dilute gases 3 . This approach has successfully predicted important characteristics of dilute granular flows 
in situations where grain-grain interactions are dissipative and thermal effects are neglig ible 

Of particular interest is the prediction of Bagnold's scaling, which states that the stress tensor scales with the 
square of the shear rate ^lj. Although kinetic theory should only hold for dilute systems, experimental ^3 an d 
numerical 0, 0, observations of Bagnold's scaling in dense systems have fostered hopes that kinetic theory can 
be applied beyond our naive expectations |16| . However, kinetic theory depends on the assumption that rigid grains, 
or "hard-spheres" , interact only through instantaneous binary collisions [171 llSj , an assumption that will eventually 
break down in dense systems. 

Recently, it has become clear ^|[23,EJ that Bagnold's scaling holds on quite general grounds, which are not limited 
by the assumptions of kinetic theory and do not require the binary collision assumption. Likewise, the structure of 
the hydrodynamic equations is largely set by conservation laws, and not by any of the assumptions that are specific 
to kinetic theory. These observations, though elementary, suggest that a direct comparison of macroscopic properties 
with the predictions of kinetic theory is not a reliable validation tool. Instead, tests of kinetic theory must rely on a 
detailed analysis of its underlying assumptions, and tailored measurements at the microscopic level. 

In this paper we use simulations to quantitatively test when the binary collision assumption is appropriate to 
describe momentum transport in hard-sphere granular flows. We begin by describing how momentum transport and 
Bagnold's scaling can be derived from macroscopic considerations, without using kinetic theory or the binary collision 
assumption. Then we present the formal derivation of momentum transport using kinetic theory. This comparison 
reveals that the microscopic formula for the stress tensor differs in each approach. We measure the ratio of the 
pressures predicted in each derivation, thereby testing a necessary condition for the validity of kinetic theory. 

I. MOMENTUM TRANSPORT WITHOUT KINETIC THEORY 

A primary objective of statistical mechanics is to provide microscopic or molecular justifications for macroscopic 
equations. At the hydrodynamic level we are interested in the macroscopic evolution of densities of the extensive 
variables mass, momentum and energy. 

In granular materials, mass and momentum are conserved, but energy is lost in each collision. An equation of 
motion for energy can still be written and involves a dissipative term. We will not focus on the energy equation in 
this discussion, but include it below for completeness. 

Before providing microscopic forms of the conservation equations, we recall the continuum forms. Mass conservation 
reads: 

where Greek superscripts represent components and repeated indices are summed. The mass density is denoted as 
p(r, t) and the fluid streaming velocity has components V a (r, t). These quantities are functions of the position r and 
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time t. Similarly, conservation of momentum reads: 

l ( ^ Q) = -^[^ y/3 + Ea/3 ]' (2) 

which defines the stress tensor H a/3 (r,t). 

Energy is not conserved in granular materials. Therefore, it is not possible to rely only on the conservation of 
energy to define hydrodynamic equations. However, it is expected that the equation for the energy can be written as 
a "conservative" part plus a dissipative part. It thus takes the form: 

ft(pe) = £;\peV° + J$ + V*V' r \+C, (3) 
where e is the energy density, Jq denotes the heat flux, and £ is the dissipative term. 



A. The Microscopic Connection 

Connecting the macroscopic fields in the hydrodynamic description to the microscopic motion of individual grains 
in a microcanonical formulation is non- trivial and has been the focus of recent studies p2l |2^ | . A microcanonical 
form of conservation equations begins with the introduction of a coarse-graining function Q, which allows us to define 
a continuous density functional: 

p(r,t) =^m i £(r-r i ), (4) 

i 

where nij and ri are the mass and position of grain i, and the sum is over all grains in the material. At the microscopic 
level, mass transport corresponds to the equation of motion for p. Taking a time derivative in Equation (gj and 
performing elementary algebra leads to 

! = -^£^« r - r (B) 

i 

where Vi is the velocity of grain i. This equation, combined with Equation Q defines the momentum density 

JV 

p(T,t)V a (T,t) = Y i m i v?g(T-T l ). (6) 

t=0 

In an analogous way, the hydrodynamic equation for momentum flux can be determined by taking the time derivative 
of the momentum density. The exact microscopic expression for the stress tensor is then determined by rewriting the 
resulting expression in terms of a divergence. This process is strictly algebraic and does not require any assumptions 
regarding the inter-particle forces It yields: 

N 1 N c i 

Z a ' 3 (T t t) = 1 £S(T- ri )m i (v?-V a )(v?-Vf , ) + - J2 <l\ / Sk-ri + sv^ds, (7) 
i=0 {i,j}=0 J ° 

where the sum is over all iV c contacts between grains, cry = ri — rj, and F°- is the contact force between grains i and 
3- 

The first term in Equation (JJJ represents how the velocity fluctuations of individual particles creates stress, and 
the second term gives the contribution from inter-grain forces. We will be most interested in the second term, which 
is often called the static contribution to the stress tensor and is denoted by 

1 Nc r 1 

Sf(r, *) = 5 £ / 0(r - r, + aa^da. (8) 

2 {«}=! J ° 

The spatial averaging of this macroscopic variable must be conducted in a different way than the spatial averaging of 
momentum density- this is due to the algebraic manipulations which converted Equation J2J i n f t he divergence of a 
tensor. 
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B. The Origin of Bagnold's Scaling 

Equations © and J7J), combined with Equation (0), comprise the exact relations for transport of momentum in 
granular materials (and any other particle or molecular materials) . The microscopic relations depend on the properties 
of every particle in the material and do not require any assumptions from kinetic theory. 

In most situations, it is impractical to determine the macroscopic variables from Equations © and J7J), since the 
precise motion of each particle is not known. In this case constitutive relations are postulated that relate the stress 
tensor to the velocity field. These constitutive relations can be tested using simulations, where the properties of each 
grain are known. Once a constitutive relation for the stress tensor is determined, the conservation equation can be 
solved for the velocity field. 

A granular simple shear flow can be characterized by the requirements that dV x /dy = 7, where 7 is the shear rate, 
and V = V x x . We limit the discussion here to properties of perfectly rigid grains, which corresponds to the limit 
where the time-scale introduced by the shear rate 7 -1 is much larger than the transport time of elastic waves. This 
hard-sphere limit has been the primary focus for kinetic theories of granular flow 0, Q . 

In the case of perfectly rigid grains, the velocity dependence of the constitutive relation for the stress tensor can be 
determined |20j |. Since the only time-scale is provided by the shear rate, the system obeys a strict invariance which 
applies to rigid grains in all flow regimes: doubling the shear rate does not alter the trajectory in configuration space, 
the system simply proceeds twice as fast. Therefore, the first term in Equation (J7J for the stress tensor must scale as 
the square of 7. In addition, upon increasing the shear rate, the inter-grain forces will scale as the square of the shear 
rate, since they are related to grain accelerations. Therefore the static stress from Equation JHJ must also depend on 
the square of the shear rate. 

This non-Newtonian behavior of the stress tensor has been observed in both experiments and simulations fl3l 
IT3, ITH l2lj of realistic granular flows, and is referred to as "Bagnold's scaling" . Simulations of granular materials in 
simplified geometries [24[ also suggest that most realistic granular flows occur in the rigid grain regime where, through 
the time invariance, Bagnold's scaling must hold. Therefore, considering the properties of perfectly rigid grains gives 
insight to the dynamics of realistic granular flows. 



II. MOMENTUM TRANSPORT WITH KINETIC THEORY 



The kinetic theory of frictionless, hard-sphere granular materials leads to a continuum hydrodynamic description, 
beginning with assumptions about the microscopic interaction between grains |2ll£| Since mass and momentum 

are conserved, the transport equations are the same as in Equations Q an d @- However, because kinetic theory is 
based on certain microscopic assumptions, the value of the stress tensor in Equation J5J) is determined within these 
assumptions. 

Recently, a formal construction of the transport equations in the hard-sphere kinetic theory framework has 
emerged jy, El El fl8| - This derivation is an extension of the kinetic theory of classical gases (or dense gases) 
where energy is conserved. The novel aspect is the introduction of a restitution coefficient a, which approximates 
the energy dissipation in individual collisions. The first s tep in the construction assumes that the only interactions 
between grains consist of instantaneous binary collisions [TtI Il8| . In each binary collision, momentum is conserved 
and energy is dissipated via the restitution coefficient a: 

(v 2 - Vi) • <7i2 = -a(V2 - vi) • <712, (9) 

where {v' 1; v 2 } denote pre-collisional velocities and {vi,v 2 } denote post-collisional velocities of the colliding grains 
indexed by 1 and 2. The normal unit vector CT12 = i^Z^i limits the dissipation to occur in the normal components 
of velocity. By enforcing Newton's equation in each collision, the dynamical rule from Equation @ determines a 
collisional impulse 

I bc = (1 + a)/xi2 [K - v 2 ) • fiia] £i2, (10) 
where /ii2 is the reduced mass. This impulse acts equally and oppositely on each grain involved in the collision. 



A. Pseudo-Liouville Equation and the BBGKY hierarchy 

After the microscopic interactions have been specified, the formal derivation of hard sphere kinetic theory relies on 
the pseudo-Liouville equation El El . This equation determines the time dependence of the iV-particle probability 
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distribution function (pdf) f^ N \ where f( N ) is defined such that f^dvjdrj is the probability that the particle 
j E [1,N] has a position rj and velocity Vj. The pseudo-Liouville equation can be derived by directly considering the 
dynamical rule in Equation © and using the assumption that only binary collisions occur. The result is 

(|+£)/W=0, (11) 

where C is a functional of f^ N \ The exact form for C is not important in this paper, and can be found in previous 

works mm. 

Once the pseudo-Liouville equation is derived, the time dependence of the s-particle pdf /( s ) can be determined by 
integration. The function /( s ) is determined directly from by the relation 



(12) 



and by integrating Equation l|ll|) over variables OiL s +i ^ v i^ r i- The pseudo-Liouville equation yields the time depen- 
dence of /( s ) as a function of an integral over C, which can be determined without further assumption |18|. 

The system of N equations for the s-particle pdfs that results from integrating the pseudo-Liouville equation is called 
the BBGKY hierarchy, and the most important of these equations for momentum transport in granular materials is 
the pdf for s = 1. This forms the basis for a derivation of transport equations and reads pj lltilisl ] 

|+<^)/ (1) (r 1 ,v lli )=T(r 1) v 1 ), (13) 

where T is an operator that describes how the one-particle pdf changes due to interactions between grains. When the 
binary collision assumption is made, this operator only depends on the two-particle pdf. The functional form is 



T(r l5 vi) = a 2 J dv 2 J daQ(a ■ g)(a ■ g) 



a 



- 2 f {2) (ri.rx -<7,v' 1)V ;,t) - (n,n +a,v 1 ,v 2 ,<)], (14) 



where g = vi — v 2 and a = aa is the continuous variable representing the displacement vector between two grains in 
contact. 

The collision operator T is often called the Enskog collision operator, since Enskog originally derived it for classical 
gases (a = 1) using a less formal approach Q. The collision operator encodes the change in the one-particle pdf due 
to collisions that produce a particle with velocity Vi (first term) and remove a particle with velocity vi (second term) . 
The step function limits the integral to pairs of colliding grains, and the factor of (a ■ g) sets the rate of collisions. 

The derivation of T explicitly relies on the binary collision assumption. If this assumption is not made, additional 
factors must be included in T, which are proportional to higher s-particle pdfs. For example, if we assume that up 
to n grains can interact simultaneously, then terms proportional to each of the s-particle pdfs, with s < n, must be 
included. Obviously, including many more terms will make the mathematics much more difficult, and values of n 
larger than two have not been considered in detail. 



B. Determining the stress tensor 

Kinetic theory furnishes Equations i|13|) and l|14fl that describe the time dependence of the one-particle distribution 
function through an evaluation of the microscopic interactions between grains. The only assumption is that the 
interactions consist solely of binary collisions. The next step is to connect the microscopic formulation to macroscopic 
transport equations, in particular momentum transport. 

Using the convention that p(r, t) — J mf^- 1 ' (r, v, t) dv, where m is the average grain mass, the momentum density 
is related to the one-particle pdf through the relation p(r, t) V(r, t) = j mvf^ 1 ' (r, v, t) dv. Therefore, by multiplying 
each side of Equation i|13[l by miVi, and then integrating over vi, we obtain an equation for the time dependence of 
the momentum density. This also results in an equation for the microscopic form of the stress tensor, as predicted by 
kinetic theory [3, C3 

Y, al3 {v,t) = J dvm(v a -V a )(v -T/ /3 )/ (1) (r,v) (15) 
+ 1±^ m CT 3 j dv^vzdaQia -g) 2 a a d f3 dsf^[r-(l-s)a,r + sa,v 1 ,v 2 }. 
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As in the more general case (Equation (JJJ), kinetic theory predicts that the stress consists of two contributions: one 
from the velocity fluctuations of individual grains and one from contacts between grains (which we earlier called the 
static stress). 

In the case of kinetic theory, the prediction for the static stress only includes contributions from binary collisions. 
As we defined static stress based on Equation J7J, we refer to the second term in Equation (|15l) as the "collisional" 
stress, and denote it as 

Kc (r, *) = J dv l( iv 2 dae(a ■ g) (a ■ g) 2 a a ^ J dsf™ [r - (1 - a)*, r + sa, Vl , v 2 ] . (16) 

The subscript "be" reminds us that this is the static stress predicted by kinetic theory, where only binary collisions 
are considered. 

Equation l|15l) also reveals that Bagnold's scaling holds for the stress tensor derived from kinetic theory. Each of 
the terms scales with two factors of a velocity, and the only velocity scale is given by the shear rate. Thus the stress 
tensor must be proportional to the square of the shear rate. This is not surprising, since any shear flow of rigid grains 
(whether or not the assumptions of kinetic theory apply) will obey Bagnold's scaling. Therefore, although kinetic 
theory is consistent with Bagnold's scaling, we must remember that the observation of Bagnold's scaling in a granular 
flow does not imply that kinetic theory holds. 



III. NUMERICAL TESTS 



The form of the momentum transport equation depends only on the conservation of momentum and therefore holds 
whether or not kinetic theory is used. However, the definition of the stress tensor is dependent on the microscopic 
assumptions made in kinetic theory. Using a kinetic theory approach we would predict Equation i|15|) . and using a 
general approach we would predict Equation @. In this section we test whether the stress tensor predicted from 
kinetic theory gives the same numerical value as that predicted by the general approach. This is a necessary condition 
for kinetic theory to make accurate predictions. 

To make the comparison we simulate a 2D polydisperse simple-shear flow of rigid grains using the Contact Dynamics 
(CD) algorithm |26l l27|. The polydispersity of the grains is characterized by an equal probability to have a grain 
diameter within the range a ± A, with A/a = 0.26. The motion of each grain is set by Newton's equations, with 
the forces between grains determined by enforcing the collision rule in Equation @ for each pair of contacting 
grains. A set timestep dt = 1CP 5 is used, and all collisions that occur in the timestep are considered simultaneously. 
Shear flow is produced using Lees-Edwards boundary conditions along with the SLLOD equations of motion. More 
details on our specific algorithm can be found in 2JJ, where we also observe that (i) Bagnold's scaling always holds, 
(ii) the polydispersity restricts crystallization, and (iii) the Lees-Edwards boundary conditions guarantee translation 
invariance. 

The CD algorithm is the natural choice for this particular numerical test. It guarantees that the grains are treated 
as perfectly rigid, so we need not worry about the value of the stiffness, as in a soft-sphere simulation [3 HE HH- 
Additionally, because the CD algorithm uses a constant timestep, it can simulate perfectly rigid granular flows at any 
density (below random close packing) and is not limited by the assumption that only binary collisions occur, as is the 
case for event driven simulations 28, 301. 



A. Comparing stress tensors 



To compare the stress tensors, we determine averages over the entire simulation cell. This simplifies the analysis 
since the averaging function is then Q = A -1 , where A is the area of the simulation cell. Additionally, because the 
simulations are translationally invariant, the positional dependence of the one and two-particle pdfs just provides a 
factor of A~ l . 

As a result, we see that the first term in Equations J7|) and l|15|l are identical. This is because the integral 
J G?v/W(r,v) = J dvA^ 1 /(v) plays the role of averaging over all grains, and can be replaced by a sum over grains. 
This reproduces the first term in Equation Q. 

Using the same simplifications, the second terms in the stress equations yield 

1 N " 
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FIG. 1: The ratio of average contact force to the average binary contribution. Values equal to one correspond to restitution 
coefficients a and packing fraction v where kinetic theory is able to make accurate predictions. 
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v af3 _ 

^bc - 



2A 



dvidv a d*0(<7 • g)(a • g )/ (2) [ Vl , v 2 ]a a l£, 



(18) 



where we have inserted the binary impulse from Equation IjlQI) into the collisional stress, with the reduced mass replaced 
by one-half the average mass. Equation (|1 T|) is the total static stress, and the Equation l|18(l is the static stress predicted 
by kinetic theory using the binary collision assumption (the collisional stress). Realizing that a 1 J dv\dv^daQ [a ■ 
g)(<7 • g)/' 2 ) is the collision rate [U, then the collisional stress can be written as 



v a/3 _ N c 



be 



2A 



(19) 



where the first term accounts for the collision rate (the number of collisions N c per simulation timestep dt) and the 
second term gives the average value of the impulse in each collision, as determined in Equation (|10|) . This can be 
transformed to a sum over collisions written as 



1 

2~4 



(Fbc) 



(20) 



where Fb c = Ibc/<^ is the total force exerted through binary collisions in a single timestep dt of the simulations. 

With the derivation of Equations H17|) and H20|) , we have succeeded in re- writing the static and collisional stress in 
forms that can be easily computed in simulations. We see that the form of the collisional stress is exactly the same as 
the static stress, with the total contact force replaced by the binary collision force. This is because the only assumption 
that kinetic theory makes to derive the stress tensor in Equation (|15fl is the binary collisions assumption. 

Therefore, in order to test the accuracy of kinetic theory, we simply need to measure the magnitude of the average 
contact force F and compare it to the average magnitude of the collisional force Fb c - A necessary test of the validity 
of kinetic theory is that this ratio is equal to one. In Figure we present data for the ratio {F)/(F\ JC ) for all values 
of the restitution coefficient and packing fraction we have investigated. The packing fraction v is equal to the area 
of grains divided by the total simulation area. The average was taken over 2000 timesteps in steady state. The 
different symbols correspond to different resitution coefficents, with the higher curves corresponding to decreasing a 
(i.e. increasing dissipation) and therefore increasing inaccuracy of predictions based on kinetic theory. 

We see from Figure ^ that for many combinations {a, v}, the ratio of forces is equal to one and kinetic theory is 
applicable. However, for larger values of v and smaller values of a, the ratio becomes larger than one and kinetic 
theory can not be used to predict the stress tensor. Notice that the failure of kinetic theory does not depend on any 
approximation scheme used to calculate constitutive relations for the stress tensor, rather it is fundamentally due to 
the assumption of binary collisions, which is made at the beginning of the formal derivation of kinetic theory. 
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FIG. 2: Contours of constant (F)/{Fb c ) — 1, with the value written beside each curve, for all restitution coefficients a and 
packing fractions v we have investigated. The different contours correspond to values of a and v where error induced by the 
binary collision assumption is 5% (0.05), 10% (0.1), 25% (0.25), 50% (0.5), and 100% (1). 



One of the more surprising conclusions of the data in Figure ^ is that the validity of kinetic theory depends on the 
pair of variables {a, v\ and not just the density. For flows at high density, a high restitution coefficient tends to keep 
grains well spaced so that only binary collisions occur; for flows at low restitution, collisions tend to cluster grains 
that have previously interacted, causing kinetic theory to fail even at low densities. 

Kinetic theory is only useful when it accurately predicts transport of momentum. In Figure |2 we plot contours of 
(F) I (Fbc) — 1 in the phase space of 1 — a and v. A plot of this quantity, which is equal to zero if only binary collisions 
occur, allows us to easily quantify the error induced by the binary collision assumption. The first thing to notice 
about the phase plot is that there is a wide range of restitution and density where the error induced by the binary 
collision assumption is less than 5%- this corresponds to the area below the curve labeled 0.05. However, depending 
on the required numerical accuracy, the binary collision assumption begins to make poor predictions at large density. 
Since most natural flows occur above a packing fraction of 0.75, the applicability of kinetic theory to natural flows is 
quite limited. 

The numerical results in Figure|3provide a quantitative test of a necessary condition for kinetic theory to be applied 
to granular flows. It is necessary that the collisional stress tensor, which is predicted by kinetic theory where only 
binary collisions are assumed to hold, be equal to the actual static stress tensor. This is only achieved for a limited 
amount of phase space. In the other regions, high dissipation and high density causes clusters of interaction grains to 
form, thereby limiting a kinetic theory approach. 

IV. CONCLUSION 

We have presented a general derivation of the momentum transfer equation and compared it to the kinetic theory 
derivation for rigid granular flows. The form of the transfer equation is based on conservation of momentum and 
therefore the same in both derivations. However, the predictions for the value of the stress tensor varies. Kinetic 
theory is accurate when the two stress tensors are equal, which occurs when the average contact force is equal to 
the average collisional force and interactions between grains are strictly binary. We have quantitatively identified the 
range where kinetic theory is invalidated by considering a ratio of the total and collisional contact forces, and we find 
that kinetic theory is unable to make accurate predictions for a certain range of density and restitution coefficient. 

In the regions where kinetic theory fails, clusters of contacting grains form that interact through force chains. The 
presence of force chains becomes the major contributor to the contact forces between grains, and these force chains 
must be incorporated in theory to properly characterize momentum transport in granular materials. 
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